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ABSTRACT 

The differentiation matrix for a Daubechies-based wavelet basis defined on an in- 
terval will be constructed. It will be shown that the differentiation matrix based on 
the currently available boundary constructions does not maintain the superconver- 
gence encountered under periodic boundary conditions. 
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1 Introduction 


Superconvergence is a property of wavelet methods, Daubechies, see [4], and spline, 
see [10], based, when periodic boundary conditions are assumed. That is, for Daubechies- 
based methods under the assumption of periodicity it was proven in [6] that the 
differentiation matrix is accurate of order 2 M, even though the highest degree poly- 
nomial which can be reconstructed exactly in the wavelet subspace is of degree M — 1, 
where M is the number of vanishing moments of the wavelet. For spline-based meth- 
ods, again assuming periodicity, it was proven in [7] that the differentiation matrix 
is accurate of order 2 n + 2, even though one can only construct splines of order n 
exactly in the underlying subspace. For Daubechies and spline-based wavelet systems 
this, roughly, doubling of the differentiation accuracy is known as superconvergence, 
a name borrowed from the finite element literature. 

When periodicity is no longer assumed and one defines wavelets on an interval 
one of the goals when building boundary functions should be to maintain the su- 
perconvergence of the differentiation matrix across the entire interval. It has been 
proven by Gottlieb, et. al. [5], that the superconvergence encountered in finite el- 
ement methods under the assumption of periodicity can not be maintained on the 
interval if there are characteristics leaving the domain. Furthermore, it was shown in 
[7] that for spline-based wavelet systems that super convergence is lost at the bound- 
aries when the boundary functions are truncated B-splines. The goal of this paper 
is to explore the accuracy of the differentiation matrix for Daubechies-based wavelet 
systems defined on an interval. In particular, the differentiation accuracy of the dif- 
ferentiation matrices constructed from the currently available boundary constructions 
for Daubechies wavelets defined on an interval will be explored. The differentiation 
matrix for the Daubechies 4 wavelet basis using the boundary functions defined in 
[3] will be calculated explicitly and the accuracy will be found. From this explicit 
example one can see a necessary condition for superconvergence to be maintained 
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with other boundary constructions. 


We begin by defining the differentiation matrix, P, as the product of three ma- 
trices: V = CDC. The three matrices are defined as follows: 

• The first matrix C is a quadrature matrix mapping from samples in the physical 
space to approximate scaling function coefficients at the finest scale: s — C f. 
C can be derived from the moments of the scaling function. (7, also, can simply 

A 

be the inverse of the third matrix: C = ((7) . 

• The second matrix D maps from the scaling function coefficients of a function 
to the scaling function coefficients of the derivative of the same function. For 
convenience 1 will, henceforth, refer to this matrix as the ‘derivative projection’ 
matrix. 

• The third matrix C is a quadrature matrix mapping from approximate scaling 

^ A 

function coefficients to approximate point values in the physical space: f — Cs. 
C can be derived from samples of the scaling function or can simply be the 
inverse of the first matrix: C = C~ l . 

This paper contains the following sections. 

1 . Introduction 

2. Boundary Functions: 

An outline of the boundary functions constructed in [3] is given. 

3. Derivative Projection Matrix D : 

The entries of D are defined and calculated. D has a banded structure due 
to the local support of the basis functions, but D is not antisymmetric. The 
number of nonzero entries in each row of D bounds the maximum differentiation 
accuracy for the corresponding row in the differentiation matrix. Furthermore, 
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the number of nonzeros entries in the first few and last few rows of D are related 
to the overlap of the boundary functions. 

4. Approximating Scaling Function Coefficients: 

An outline of scaling function coefficient approximation is given. 

5. Quadrature Matrix C Using Moments: 

C performs the mapping, C : / — ► s. The D\ wavelet has two vanishing 
moments which yields well-conditioned 2-point quadrature matrices and ill- 
conditioned 3-point quadrature matrices. 

6. Quadrature Matrix C Using Samples: 

A A ^ A 

C performs the mapping, C : s — ► /. For the D 4 wavelet the matrix C has 
three entries in each row and is banded. C is, however, ill-conditioned. 

7. Examples and Order of Accuracy: 

The matrix D is fixed. There are, however, many choices for the quadrature 
matrix producing many versions of T > . In all cases superconvergence is lost at 
the boundaries. Furthermore, choosing among the quadrature matrices which 
not ill-conditioned one gets a differentiation matrix which is essentially full. A 
banded V can be constructed from a combination of the two types of quadrature 
matrices C and (7. The bandwidth from such a construction is large, 9, for 
a D4 wavelet and, again, the superconvergence is lost at the boundaries. In 
addition, all of these matrices display directional dependency. This directional 
dependency can be ‘fixed’ by reversing the coordinate system and averaging the 
appropriate differentiation matrices. 

8. Conclusion: 

The differentiation matrices for wavelets on an interval using the boundary func- 
tions defined in [3] do not maintain the superconvergence encountered when 
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periodic boundary conditions are assumed. By simply counting the overlap of 
boundary functions in other boundary constructions it can be seen that all cur- 
rently available boundary constructions fail to maintain the superconvergence 
encountered under periodic boundary conditions. 
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2 Construction of Boundary Functions 

This section is a brief outline of the boundary function construction proposed by 
Cohen, Daubechies, and Vial, see [4]. 

The goal is to build a wavelet basis on an interval where the scaling functions and 
wavelets away from the boundaries are the usual Daubechies scaling functions and 
wavelets. At the boundaries, boundary scaling functions are constructed such that 
polynomials of order up to the number of vanishing moments of the wavelet can be 
reproduced exactly across the entire interval. 

The boundary function construction begins by building independent, but not or- 
thogonal, functions as follows: 

2N-2 / \ 

#*(*)= £ (l U(* + n-W+l), (1) 

where <f>(x) is the usual Daubechies scaling function and N is the number of van- 
ishing moments of the associated wavelet. Note that these functions are compactly 
supported, and their supports are staggered, i.e., 

$upp((f> k ) = [0,2iV — 1 — A:]. (2) 

The staggered support yields independence, and the boundary functions are then 
defined by simply orthonormalizing these functions by a Gram-Schmidt method. 

Coefficients h L and h R can be found, see [4], such that the boundary functions are 
defined recursively as follows beginning with the left-hand side boundary functions: 

N - 1 N+2k 

<p L k {x) = v/2 £ ti tl tf(2x) + V2 £ h km <f>(2x - m), (3) 

/=0 m—N 

and the equation defining the right-hand side boundary functions is, 

N-l N+2k 

K(x) = \/2 £ (2x) + V2 £ + m + 1 ). (4) 

/=0 m — N 

The following figure is a plot of the four boundary scaling functions for the 


Daubechies 4 wavelet. 




Figure 1 : The boundary functions for the Daubechies 4 scaling function 
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3 The Derivative Projection Matrix D 

Let us recall the origin of the matrix D. First, we begin with a function f(x) 6 
L 2 (R). Next, we want to approximate, or discretize, this function in a scaling function 
subspace Vo C L 2 (R), 

p v 0 f(x) = s k b k (x), (5) 

k = 0 

where {bk(x)} denotes the bases functions of this finest scale discretization and d 
is the dimension or number of degrees-of-freedom of the subspace. In the case of 
wavelets with periodic boundary conditions bk(x) — <j>k(x ), the usual Daubechies 
scaling functions, for all k £ — 1}. In the current scenario, however, of 

wavelets on an interval b^(x) is different at the boundaries: 

b k(x) = <l>k( x ) 


for k = 0 , N — 1, 


bk{x ) = 4>k{x) 


for N < k < d — 2 N + 1 and, 

b k( x ) = 4>k( x ) 

for k = d — 2 N + 2 ,...d — 1, where N is the number of vanishing moments of the 
wavelet. Under the orthonormality of the current basis set, < 6jt, 6/ >= 6k, i, where, 


<g,h>= 


J g(x)h(x)dx , 


and denotes integration over the interval, we get, 


( 6 ) 


Sk =< &*,/>• 


( 7 ) 


We now differentiate the approximation Pv 0 f(x) to get, 


d_ 

dx 


Pv Q f(x) = Yl S k b k(x), 


k = 0 


( 8 ) 
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where b(x) = -^b(x). The derivative takes Pv 0 f(x) out of Vo. Projecting back into Vo 
we get, 

Pv ‘Ji PvJ{z) = S < ^ Pv ° f ’ b ' > 6 'W- < 9 ) 

or, 

Pv oT PvJ(x) = < b k ,b t > bi(x). (10) 

ax 1=0 k =0 

/ 

The elements < 6*, 6/ > comprise the derivative projection matrix D . If we let s 
denote the vector of the scaling function coefficients s*, for k = 0 , — 1, then D 

maps from the scaling function coefficients of a function at the finest scale to the 
scaling function coefficients of the derivative of the same function: 

D : s s. 

D , in effect, contains all the derivative information in a wavelet expansion. The 
elements of D contain the numerical values of the derivative of each scaling function 
projected onto all scaling functions in the basis. Due to the compact support of 
the bases functions, the matrix D has a band diagonal form. The following matrix 
illustrates the form of D for the D\ wavelet basis with boundary functions defined as 
in the previous section: 


Po,o 

Pi, 0 

a 2,0 

0 

0 

0 

0 

0 

0 

0 ■ 

P 0,1 

Pi, 1 

«2,1 

Q 3,l 

0 

0 

0 

0 

0 

0 

00,2 

01,2 

r 0 

n 

r 2 

0 

0 

0 

0 

0 

0 

01,3 

r - 1 

r 0 

n 

r 2 

0 

0 

0 

0 

0 

0 

r- 2 

r - 1 

r 0 

ri 

r 2 

0 

0 

0 

0 

0 

0 

r- 2 

r - 1 

r o 


r 2 

0 

0 

0 

0 

0 

0 

r- 2 

r - 1 

r o 

r i 

0?A 

0 

0 

0 

O' 

0 

0 

r- 2 

r-\ 

TO 

01,5 

0£s 

0 

0 

0 

0 

0 

0 

n R 

a 4,l 

<1 

Phi 

Po, 1 

0 

0 

0 

0 

0 

0 

0 

<*S,0 

Pi, 0 

Pofi . 


where the coefficients r are the derivative projections due to scaling function inter- 
action, the coefficients p are the projections of the derivatives of boundary functions 
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projected onto boundary functions, and the coefficients a and f3 contain the interac- 
tions between derivatives of boundary functions and scaling functions. 

In the following subsections the entries of D will be found. 

3.1 The Equations for the Derivative Projections 

To find the elements of D one needs to find the interaction of the derivative of each 
basis function with every other basis function. Due to the local support of the bases 
functions, the left-hand side boundary functions and the right-hand boundary func- 
tions due not interact. (Assuming, of course, that one does not choose an unreason- 
ably small domain relative to the support of the scaling functions.) 

Each entry in the matrix D must be found separately, but the calculations are 
straightforward. A sampling of the relevant derivations will be given below. 

3.1.1 Derivation of the LHS Equations 

Recall that the LHS boundary functions <j>k(x) are generated by, 

N- 1 N+2k 

<f>k{x) = y/2j2 (2*) + ^2 E rt tm <f>(2x - m). (11) 

/=0 m=N 

Differentiate the above equation to get, 

t /V-l JV+2* 

i-1 4(as) = 2^2 £ +2^2 £ hj;„<H2x - m). (12) 

dx 1=0 m=N 

Multiply <j>k( x ) by <j>y{x) and integrate to get, 

f <j>%(x)<f>p(x)dx = 4 E h k,t h i,i fj>j'{2x)<f>f(2x)dx (13) 

Jl 1=0 i=o JI 

N+2kN-l . 

+ 4 E E h i,i h k,rn \. ti (2x)<f>{2x - m)dx 

m=N «'= 0 Jl 

N- 1 N+2 p 

+ 4 E E h k,l h ™ L - q)dx 

1=0 q=N JI 

N+2kN+2p . 

+ 4 E E h k,m h v, q I. <K 2 * ~ m)<f>(2x - q)dx. 

m=N q=N JI 
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Rename the integrals as, 


pir = l <t>i(x)<t>i(x)dx, 

*i,i = J j ^>{x-m)<j>f(x)dx, 
fa = J j <l>i(x)<f>{x-q)dx, 
r m ,q = f^{x - q)dx, 


(14) 

(15) 

(16) 
(17) 


to get the following system of equations: 

N- 1 TV- 1 


TV-1 AT+2 k 

l l 2 pir = EE4W + E 

/ i i=0 m=N 

N-\ N+2 V N+2k N+2p 

E E W.A, + E E '■tAw 

1=0 q=N m=N q=N 


(18) 


The numerical values of the coefficients {r} were found in [2], and the numerical 
values of the coefficients p L , a L , and /? L will be found in the following subsections. 

3.2 The Coefficients for i = 0, 1 

In this subsection we will find the corner diagonal elements of D. These elements 
represent the projection of the derivative of each boundary function onto itself. 

3.2.1 The Coefficients p% 0 and pf, 

For simplicity, allow the lower limit to be 0 and the upper limit be oo in the definition 
of Pk,k- 

too , 

Pk,k = / <t>i{x)<t>i{x)dx. 

Jo 

Apply integration by parts to this equation to get, 


Pk,k = (^( I )) 2 |S° - / <f>k(x)<f>k(x)dx, 
Jo 


(19) 


or. 


2 pU = (^(^)) 2 !~- 


( 20 ) 
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( 21 ) 


<f> k (x) is 0 at x = oo leaving only the value at the lower limit, 

Pk,k ~ ~ 2 

That is, p kk is determined by the value of the boundary function at the boundary. 
Before evaluating the boundary functions at the boundary it will be shown that pfc 0 
and p\ x are related by a constant multiple. 

3.2.2 The Relationship Between p^ Q and p\ x 

In equation (11) set x = 0 to get, 


1 2+2 it 


4>k( 0) = V2 1 >^( 0 ) + V2 Y, h \ 

/=0 m= 2 

(22) 

For k = 0 we get, 


4 ( 0 ) = ^ 2 « o 4 ( 0 ) + K. , 4 («)). 

(23) 

which yields, 


<f>i( 0 ) = c L <f>o(0), 

(24) 

where 

L 1 — \/2/lo,o 

y/2 

(25) 


Recall from above that the p ktk 's can be found from the boundary functions evaluated 
at the boundary. Combine this information with equation (24) to get, 


A, = (c'OVoV < 26 ) 

That is, we only need to evaluate one of the boundary functions at the boundary 
in order to know pQ 0 and p\ x . The following subsection will evaluate the boundary 
function (f> q(x) at x = 0. 

3.2.3 Evaluation of the Boundary Function <f>Q (x) at x = 0. 

In this subsection the straightforward evaluation of the boundary function <fio(x) at 
x = 0 will be given. 
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In equation (11) let k = 1 and x = 2 to get, 


tf(2) = v/2(A^(l) + k^O))- (27) 

The values of the scaling function (f>(x) at the above integers are <^(0) = 1/2(1 + \/3) 
and <f>{\) = 1/2(1 — \/3), see [9]. Using these numbers and the coefficients provided in 
[3] we get <f>}( 2) = —.3654971039. Similar to above, let k = 0 and x = 1 in equation 
(11) to get 

*o0) = V2h^(2) + y/2h^{ 0), (28) 

and evaluate to get = —1-1265992786. In this way, one can get the following 

values for the boundary functions: 1) = 1.4058555636, 1/2) = .4123639530. 

Now we have two values for 4>q(x) at x = 1/2 and x = 1. Furthermore, the two 
boundary functions are lines between x = 0 and x — 1 . We can, therefore, find 
the value of <f> q(x) at x = 0: <j> q( 0) = 1.9513271847. As outlined in the previous 
subsection we now have sufficient information to find Pq 0 and p\ x : 

Po,o = -1-9038388908, 

and 

pfj = -.04295212263. 

The calculation of the remaining elements of the derivative projection matrix is 
straightforward and will not be given. The next subsection provides the complete 
matrix D. 

3.3 The Complete Matrix D 

This section of the paper will be concluded by giving a 10 by 10 version of the matrix 
D. Only four digit accuracy will be given for the corner entries p, a, and (i. The 
entries r are known exactly. 


D = 


12 



1.9038 

.9444 

-.2565 

0 

0 

0 

0 

0 

0 

0 

1.5163 

-.0430 

.6752 

-.0832 

0 

0 

0 

0 

0 

0 

2565 

-.6752 

0 

2 

3 

__L 

12 

0 

0 

0 

0 

0 

0 

.0832 

2 

3 

0 

2 

3 

1 

12 

0 

0 

0 

0 

0 

0 

j_ 

12 

2 

3 

0 

2 

3 

12 

0 

0 

0 

0 

0 

0 

1 

12 

2 

3 

0 

2 

3 

i 

12 

0 

0 

0 

0 

0 

0 

J_ 

12 

2 

3 

0 

2 

3 

-.0765 

0 

0 

0 

0 

0 

0 

J_ 

12 

2 

3 

0 

.5825 

-.0397 

0 

0 

0 

0 

0 

0 

.0765 

-.5825 

.0899 

.3150 

0 

0 

0 

0 

0 

0 

0 

.0397 

-.7936 

.6369 


Note that, as expected, away from the boundaries D is the same as the derivative 


projection matrix which was found for periodic boundary conditions, see [6]. Also, 
note D is not antisymmetric, but is banded. D is not in itself the differentiation 
matrix, but D is the component of the differentiation matrix, T) — C 1 DC , which 
does the differencing. Later in this paper it will become obvious how one can simply 
‘look’ at D in order to tell the maximum differentiation accuracy at the boundaries. 
The next section begins the discussion of choosing a quadrature matrix, C . 
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4 Approximating Scaling Function Coefficients 

Once again, let us work with only the D 4 wavelet. Recall, as constructed in [3], 
there are two boundary functions at each end of the interval, as well as, the usual D 4 
scaling functions away from the boundaries. Let N denote the number of degrees-of- 
freedom of the approximation space, Vo- Our function f(x ) € L 2 (R) is then discretized 
as follows (In order to simplify the presentation I will let the ‘interval’ be [0, N], 
and there will be no overlap between left-hand boundary functions and right-hand 
boundary functions.) 

First, the left-hand side, 

sf = f <f>i(x)f(x)dx , (29) 

Jo 

then the right-hand side, 

= / <f>^(x)f{x)dx, (30) 

Jo 

where i = 0, 1. Next, discretize in the middle with the usual scaling functions, 

r N 

Sk = I Mx — k)f(x)dx, (31) 

J 0 

for k — 2, ..., N — 4. 

Now we will assume that f(x) is defined on a discrete evenly-spaced grid, /(x<) 
for i = 1 We must, therefore, approximate the above three integrals by an 

appropriate quadrature matrix. 

In the following two sections two types of quadrature matrices will be consid- 
ered. The first type is constructed using the moments of the scaling function and is 
represented by, 

s = Cf. (32) 

It will be seen that two-point quadrature matrices for this method are well-conditioned, 
whereas three-point quadrature matrices are ill-conditioned. 
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The second type of quadrature matrix is constructed from samples of the scaling 
function and boundary functions and is represented by, 

/= Cs. ( 33 ) 

It will be seen that this type of quadrature matrix is always ill-conditioned. 
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5 Quadrature Matrix C Using Moments 


The quadrature matrix C maps from evenly-spaced point values, /(x,-), of a function 
to the approximate scaling function coefficients at the finest scale, s. The quadrature 
matrix depends on three quantities: the number of vanishing moments of the wavelet, 
the grid that the function f(x) is defined on, and which grid points are used to 
approximate each scaling function coefficient. 

5.1 Moments of Scaling and Boundary Functions 

The moments of the D 4 scaling functions were calculated in appendix A of [6]. The 
first 3 moments are, Mo — T Afi = .6339745962, and M 2 = .4019237886. 

We will concentrate here on calculating the moments of the boundary functions. 


5.1.1 Moment 0 of Left-Hand Side Boundary Functions 


The 0-th moment of the boundary function k is, 

- q© 1 j»oO 2+2 k j»oo 

/ ti(z)dx = y/\ 2j2 h tt ^(2x) + V 2 £ h L k<m / <f>(2x-m). (34) 

J 0 l =0 Jo m—2 J ° 

Now, adopt the notation m° k = f cf> k (x)dx and let y = 2x where appropriate in the 
above integrals to get, 


1 2+2 k m oo 

V2m° k = + £ h L km j <f>{x)dx. (35) 

/=0 m—2 J ~ m 

The integral f™ m <f>(x)dx is 1 since m > 2, and the scaling function, as defined here, 
is supported on [—1,2]. Equation (35) defines the following 2 by 2 system: 

(36) 


sr- 

pt* 

0 

1 

5, 


m 0 



~ho2 j 

L 


m j 


h\ t 2 

^i-T 

1 

CO 

1 


The solution to the above system is 

mg = .3620520589, 
m? = 1.001445402. 

In a similar fashion all moments can be found for the left and right hand side 
boundary functions. 
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5.2 The Choice of Grid Points 


In this subsection the affect of the choice of different grid points on the quadrature 
matrix will be illustrated. For simplicity, in this subsection we will work with only 
8 bases functions and 8 grid points. The interval will be [0,8]. When choosing our 
grid points to approximate the scaling function coefficient s, we must use only the 
grid points which lie within the support of the basis function 6,(x). That is, the inner 
product that is being approximated is, 
fS 

Si ~ 

Jo 

and the data that is available to approximate s, must be limited to function values 
defined within the support of 6,(:r). The support of each of the eight bases functions 
is, 

supp{(f)\ d } = [0,2], 
supp{4> b 2 d ] = [0,3], 
supp{<f> 3 } = [1,4], 
supp{<j> 4 } = [2,5], 
supp{<f) 5 } = [3,6], 
supple) = [4,7], 
supp{<t> bd ) = [5,8], 
supp{<t> 8 d } = [6,8], 

where the superscript ‘M’ denotes a boundary function. 

The following subsections will consider a few of the possible quadrature matrices 


17 



corresponding to the following two grids: 


and 


' .5 ‘ 

1.5 

2.5 

3.5 

4.5 

5.5 

6.5 

7.5 


y = 


o 

8/7 

2 * 8/7 
3*8/7 
4*8/7 
5*8/7 
6*8/7 
8 


5.2.1 The Entries of the Quadrature Matrices 


( 37 ) 


(38) 


Now, one must use the scaling function and boundary function moments and at least 
two grid points contained within the support of each basis function in order to build 
the quadrature matrix. The first row of the quadrature matrix will estimate the basis 
function coefficient for the boundary function with the smallest support, <f>\ d . Suppose 
we have chosen grid x. We know that for the D 4 wavelet that constants and lines can 
be represented exactly in the wavelet subspace Vo- There are only two grid points 
within the support of <f>\ d . Therefore, the solution to the following linear system will 
yield the two coefficients c^i and C 12 in the first row of the quadrature matrix. The 
linear system is, 


fi <f>\ d {x)dx 


' 1 1 

c i,i 

f] (f>\ d (x)xdx 


.5 1.5 

c \,2 


In a similar manner all the entries of a quadrature matrix can be found once one 
has calculated the moments of the basis functions and chosen the grid points. The 
following subsection will illustrate a few of the quadrature matrices associated with 
various grid choices. 
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5.2.2 All Grid Points Contained in Scaling function Support 


As a first possible choice of grid points, let us use every grid point that is contained 
within the support of the respective basis function. Such a quadrature matrix will 


the following structure: 



Cl.l 

Cl ,2 

0 

0 

0 

0 

0 

0 

C2 , 1 

c 2,2 

c 2,3 

0 

0 

0 

0 

0 

0 

C3,2 

C3,3 

C3,4 

0 

0 

0 

0 

0 

0 

c 4,3 

C4,4 

C4,5 

0 

0 

0 

0 

0 

0 

C5,4 

C5,5 

C5,6 

0 

0 

0 

0 

0 

0 

C6,5 

Q>,6 

C6,7 

0 

0 

0 

0 

0 

0 

C 7,6 

c 7,7 

C7,8 

0 

0 

0 

0 

0 

0 

C8,7 

C8,8 


(40) 


For both of the grids x and y the matrix C^ pt is ill-conditioned. The condition number 
is proportional to N 2 , where N is the number of grid points, and inversion of 


unacceptably corrupts the differentiation matrix. 


5,2.3 2-Point Quadrature at Boundary and 3-Point Otherwise 


The quadrature matrix has the following structure: 




c l,l 

Cl ,2 

0 

0 

0 

0 

0 

0 ' 




0 

c 2,2 

c 2,3 

0 

0 

0 

0 

0 




0 

C3,2 

C3,3 

C3,4 

0 

0 

0 

0 


£r3pt 

= 

0 

0 

0 

0 

C4,3 

0 

C4,4 

C5,4 

C4,5 

C5,S 

0 

C5,6 

0 

0 

0 

0 

(41) 



0 

0 

0 

0 

C6,5 

C6,6 

C6,7 

0 




0 

0 

0 

0 

0 

C7,6 

C7,7 

0 




0 

0 

0 

0 

0 

0 

C8,7 

C8,8 


This matrix is, also, 

ill-conditioned for 

both 

grids 

x and y. 

Again, the condition 


number if proportional to N 2 . 
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A First Option Using a 

2-Point Quadrature 



C\,\ 

Cl ,2 

0 

0 

0 

0 

0 

0 


0 

C2,2 

C2,3 

0 

0 

0 

0 

0 


0 

0 

C3,3 

C3,4 

0 

0 

0 

0 

{j2 pt 

0 

0 

0 

C4,4 

C4,5 

0 

0 

0 

0 

0 

0 

0 

C5,5 

C5,6 

0 

0 


0 

0 

0 

0 

0 

C6,6 

C6,7 

0 


0 

0 

0 

0 

0 

C7,6 

C7,7 

0 


0 

0 

0 

0 

0 

0 

c 8,7 

C8,8 


This quadrature matrix is well-conditioned for both grids x and y. 


( 42 ) 


Second Option 

. for 2-Point Quadrature 

r c,,i ci ,2 oooo 

0 

0 


c 2,l 

C2,2 

0 

0 

0 

0 

0 

0 


0 

0 

C3,3 

C3,4 

0 

0 

0 

0 

c?' = 

0 

0 

0 

C4,4 

^4,5 

0 

0 

0 

0 

0 

0 

0 

c 5,5 

C5,6 

0 

0 


0 

0 

0 

0 

0 

C6,6 

C6,7 

0 


0 

0 

0 

0 

0 

0 

c 7,7 

C7,8 


0 

0 

0 

0 

0 

0 

C8,7 

C8,8 


This quadrature matrix is, also, well- conditioned for grids x and y. 


(43) 


5.3 Conclusion 

In this section a variety of structures for quadrature matrices have been shown where 
the matrices are found by using the moments of the scaling function. In general, the 
pattern is that 3-point quadrature matrices are ill-conditioned and 2-point quadrature 
matrices are well-conditioned. 

The next section will explore quadrature matrices constructed from the samples 
of the scaling functions. 
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6 Quadrature Matrix C Using Samples 


Let g(x) denote the projection of f(x) 6 L 2 (R) in Vo: 

g(x) = P v J(x) = s kh{x). (44) 

k = 1 


The quadrature formula based on samples of the basis functions can be found 

from, 

N 

(45) 


g(x,) = PvJ(xi) = Y, Skh(xi). 

k=\ 


In matrix form the quadrature formula is, 


9= Cs, 


(46) 


or let a be the approximation to s to get, 


f = Ca, 


(47) 


where for an 8 by 8 case the quadrature matrix is, 


c = 


<f>o(x i) 

^0(^2) 

0 

0 

0 

0 

0 

0 


<t>\{x 1 ) 

<t>i{x 2 ) 

</>l( X 3) 

0 

0 

0 

0 

0 


0 

Mx?) 

fa{x 3 ) 

<^2(^4) 

0 

0 

0 

0 


0 

0 

h(x 3 ) 

fa(x 4 ) 

<t>: 3(^5) 

0 

0 

0 


0 

0 

0 

<M* 4 ) 

^4(^5) 

M x e) 

0 

0 


0 

0 

0 

0 

<t>s{Xh) 
<f>s(x t 3) 
<j>s(x 7 ) 

0 


0 

0 

0 

0 

0 

^*(*«) 
<t>l{x 7 ) 
<i> 4 {xs) 


0 

0 

0 

0 

0 

0 

<t> 0 (^ 7 ) 

<t>o( x s) 


(48) 


Now C will be 


found for the two previously defined grids x and y. 


6*1 A Quadrature Matrix based on Grid x 

Finding the explicit form of the above defined C for the grid x is relatively simple since 
we only need the numerical values of the scaling functions and boundary functions 
halfway between the integers. 
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6.1.1 

The Numerical Values of the Scaling Function 


Using the definition of the scaling function supported on [0,3], 



3 

<f)(x) = \/2 5Z 2® - &)> 

k= 0 


we get 

^(5/2) = V2h 3 <f>{2) = .06698729811, 

(49) 


<f>{ 3/2) = V2{h\(j>{2) + /*2<Hl))> 

(50) 

and 

0(1/2) = v/2M(l) = -9330127019. 

(51) 


The numerical values of <j) at the integers are provided in [9]: 

<f>( 1) = 1/2(1 +V3), 

and 

<j>{2) = 1 / 2(1 - VS). 

6.1.2 The Numerical Values of the LHS Boundary Functions 

Given the numerical values of the scaling function <f>(x) on grid x one can now find 
the numerical values of the boundary functions on grid x from the expression, 

1 2+2Jfc 

<f>k(x) = \/2 E ^ E h k,mtt 2x - m )• ( 52 ) 

(=0 m = 2 

In a similar manner one can find the numerical values for the RHS boundary 
functions to get the quadrature matrix of the next subsection. 

6.1.3 The Complete Quadrature Matrix for Grid x 

An 8 by 8 example of the quadrature matrix for grid x is, 
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.4124 

.8495 

0 

0 

0 

0 

0 

0 

.2062 

-.0077 

.9330 

0 

0 

0 

0 

0 

0 

.0669 

0 

.9330 

0 

0 

0 

0 

0 

0 

.0670 

0 

.9330 

0 

0 

0 

0 

0 

0 

.0670 

0 

.9330 

0 

0 

0 

0 

0 

0 

.0670 

0 

.8561 

0 

0 

0 

0 

0 

0 

.0670 

.3270 

.4451 

0 

0 

0 

0 

0 

0 

-.1406 

.8902 


( 53 ) 


This matrix is ill-conditioned just as the 3-point quadrature formulas were for the 
quadrature matrix derived using moments. 


6*2 A Quadrature Matrix based on Grid y 

Using the definition of the scaling function, 

_ 3 

<f>(x) = V2 h k <f>( 2x - k ), 

Jt=0 

one can set up a matrix such that the eigenvector for the eigenvalue A = 1 gives a scalar 
multiple of the scaling function 4>{x) at x = if 7 for i — 1,...,20. The eigenspace for 
the eigenvalue A = 1 is, however, 2-dimensional. To avoid this complication one can 
approximate the numerical values of 4>(x) at x = i/7 for i = 1,...,20 by generating 
<f>{x) at a large number of diadic points, in this case 2 14 points were chosen. This 
approximation is sufficient to understand the character of the quadrature matrix, 
specifically if it is well-conditioned or not. Once the numerical values of the scaling 
function are known then we can find the numerical values of the boundary functions. 


6.2.1 The Numerical Values of the LHS Boundary Functions 


As before, the numerical values of the LHS boundary functions on grid y can be found 
from, 


2+2fc 


<t> L k {x) = \Z2^h^l'(2x) + V2 22 h k,m4>( 2x - m )- 


/=0 


m= 2 
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The numerical values in terms of the scaling function values on grid y are, 

(8/7) = V2(h^ #(16/7) + fcfc#2/7)) 

#(8/7) = #2(#i#(16/7) + # 2 # 2/7) + # >3 #-5/7)) (54) 

#(16/7) = v^(###l/7) + # 4 #4/7)). 

As with grid x the quadrature matrix produced using grid y is ill-conditioned. Re- 
call from the previous section that all 3-point quadrature matrices were ill-conditioned. 
I, therefore, conclude this section with the observation, not a proof, that quadrature 
matrices for wavelet bases on an interval using the construction in [3] should be con- 
structed using moments and that the number of grid points used to approximate each 
scaling function coefficient should be equal to the number of vanishing moments of 
the corresponding wavelet. 
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7 Order of Accuracy 

In this section the order of accuracy for a few of the differentiation matrices which can 
be constructed from the derivative projection matrix D and the quadrature matrices 
constructed in the previous three sections will be found. We begin with the differ- 
entiation matrices constructed using grid x followed by the differentiation matrices 
constructed using grid y. 

7.1 Accuracy Using Grid x 

Recall that grid x is comprised of the points halfway between the integers and, con- 
sequently, the spacing between the grid points is the same as the translation distance 
for the scaling functions on the finest scale. Furthermore, this grid does not provide a 
grid point at the boundary itself. This lack of a grid point at the boundary makes this 
grid an unlikely choice, but I have chosen to study the accuracy as an experiment. 
Let us begin by giving a few explicit examples of differentiation matrices. 

7.1.1 Differentiation Matrix DC\ vt , Grid x 

Note the unusual form of the following matrix. It is essentially upper triangular in 
form. This form is due to the inversion o f C^ pt . From Table (7.1.3) it can be seen that 
this differentiation matrix does not maintain the superconvergence at the boundaries. 
That is, the differentiation is lst-order accurate at the boundaries, whereas it is 4th- 
order accurate at the middle 50% of the grid points. The total differentiation accuracy 
across the entire domain is, therefore, 2nd-order. 

D = 
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' - 25.33 

42.38 

- 20.77 

3.73 

1 

1 

0 

0 

— 7 e " 4 

7 e " 3 

— 9 e -3 

2 e -5 

3 e " 3 

- 5.71 

- 1.59 

8.35 

- 1.04 

- le -3 

0 

0 

— 2 e -3 

.02 

-.02 

3 e -5 

6 e - 3 

1.05 

- 8.11 

.05 

8 

-1 

0 

0 

— 3 e -3 

.03 

-.04 

7 e -5 

.01 

0 

1 

-8 

0 

8 

-1 

0 

— 7 e -3 

.07 

-.09 

2 e ~ 4 

.03 

0 

0 

1 

-8 

0 

8 

-1 

-.02 

.16 

-.20 

3 e -4 

.06 

0 

0 

0 

1 

-8 

0 

8 

- 1.03 

-.34 

-.43 

7 e " 4 

.13 

0 

0 

0 

0 

1 

-8 

0 

7.93 

-.28 

-.94 

2 e ~ 3 

.29 

0 

0 

0 

0 

0 

1 

-8 

-.16 

9.56 

- 3.02 

3 e - 3 

.62 

0 

0 

0 

0 

0 

0 

1 

- 8.35 

3.36 

3.65 

-.99 

1.33 

0 

0 

0 

0 

0 

0 

0 

.25 

-.76 

- 8.37 

6.02 

2.86 

0 

0 

0 

0 

0 

0 

0 

- 1.62 

16.61 

- 33.66 

11.96 

6.71 

0 

0 

0 

0 

0 

0 

0 

.53 

- 5.47 

2.78 

- 3.27 

5.43 


7.1.2 Differentiation Matrix ( C % pt ) l DC% pt , Grid x 

Now we use the quadrature matrix C\ pt . But, the results are similar to the previ- 
ous matrix. The structure is unacceptable and the superconvergence is lost at the 
boundaries as can be seen from Table (7.1.3) 


V = 


- 12.21 

7.44 

10.48 

- 6.44 

.74 

0 

0 

0 

0 

0 

0 

0 

- 10.99 

- 11.15 

44.98 

- 24.53 

1.70 

0 

0 

0 

0 

0 

0 

0 

- 2.50 

- 1.01 

- 3.50 

8 

-1 

0 

0 

0 

0 

.01 

-.01 

0 

.46 

.07 

- 7.54 

0 

8 

-1 

0 

0 

0 

.01 

-.02 

.01 

0 

0 

1 

-8 

0 

8 

-1 

0 

-.03 

.02 

-.04 

.02 

0 

0 

0 

1 

-8 

0 

8 

-1 

-.01 

.05 

-.09 

.04 

0 

0 

0 

0 

1 

-8 

0 

8 

- 1.01 

.11 

-.19 

.09 

0 

0 

0 

0 

0 

1 

-8 

0 

7.97 

-.75 

-.41 

.19 

0 

0 

0 

0 

0 

0 

1 

-8 

-.06 

8.53 

- 1.89 

.04 

0 

0 

0 

0 

0 

0 

0 

1 

- 8.13 

1.14 

6.09 

-.11 

0 

0 

0 

0 

0 

0 

0 

0 

.73 

- 5.53 

- 3.12 

7.93 

0 

0 

0 

0 

0 

0 

0 

0 

-.24 

2.74 

- 16.75 

14.26 


7.1.3 Conclusion for Grid x 


Recall that in the construction of the differentiation matrix that the middle matrix D 
is fixed. The Table (7.1.3) contains the accuracy for two well-conditioned quadrature 
matrices C^ pt and C % pt . In addition, as an experiment I have combined two ill- 
conditioned matrices: one mapping from f to s and the other mapping from s to 
/• 
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Differentiation 

Grid 

Boundary 

Error 

Inner 

Error 

Total 

Error 

Matrix 

Size 

l\ Error 

Ratio 

l\ Error 

Ratio 

/1 Error 

Ratio 

{C'?')-' DC? 1 

16 

.0066 


2.86e -4 


.0024 


Grid x 

32 

.0032 

2.06 

3.20e -6 

89 

5.73e -4 

4.19 


64 

.0016 

2.00 

2.10e -9 

1524 

1.41e -4 

4.06 


128 

7.85e" 4 

2.04 

2.80e -11 

75 

3.50e -5 

4.03 


256 

3.81e -4 

2.06 

1.75e -12 

16.0 

8.53e -6 

4.10 


512 

1.69e -4 

2.25 

1.07e -13 

16.4 

1.94e -6 

4.40 

(C?*)-' DC?' 

16 

.0091 


2.71e -5 


.0026 


Grid x 

32 

.0047 

1.94 

3.13e -7 

86.58 

6.78e -4 

3.83 


64 

.0024 

1.96 

5.92e -10 

529 

1.72e -4 

3.94 


128 

.0012 

2.00 

2.80e -11 

21.1 

4.33e" 5 

3.97 


256 

6.08e -4 

1.97 

1.75e -12 

16.0 

l.lOe -5 

3.94 


512 

3.22e -4 

1.89 

1.07e -13 

16.4 

2.89e- 6 

3.81 

C s DC'i pt 

16 

.0038 






Grid x 

32 

.0019 

2.00 






64 

9.23e -4 

2.06 

3.95e -6 

4.05 

7.87e -5 

3.99 


128 

4.64e~ 4 

1.99 

9.83e -7 

4.02 

1.97e -5 

3.99 


256 

2.32e -4 

2.00 

2.45e -7 

4.01 

4.92e -6 

4.00 


512 

1 . 1 8e — 4 

1.97 

6.12e -8 

4.00 

1.22e~ 6 

4.04 


Table 1: Differentiation accuracy for Daubechies wavelets defined on an interval using 
grid points halfway between the integers. 


The noteworthy points are, first of all, that in all three cases the accuracy at the 
boundary is first order. Recall that when periodic boundary conditions are imposed 
that the differentiation accuracy is 4th order for the D 4 wavelet. This loss of su- 
perconvergence is a serious problem. Furthermore, note from the table that the two 
well-conditioned quadrature matrices maintain the superconvergence away from the 
boundary. But, in all three cases the total error of the differentiation matrix across 
the entire domain is one higher than the boundary accuracy, i.e., we get 2nd order 
differentiation accuracy. 

Note: Inner accuracy is calculated at the middle 50% of the grid points, and 
boundary accuracy is calculated at the two outermost grid points at each end of the 
interval. 
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7.2 Accuracy Using Grid y 

We now study the differentiation accuracy where we have included a boundary point 
at each end of the interval in the grid and the remainder of the grid points are 
required to be evenly-spaced across the interval. This choice of grid is a much more 
likely choice than grid x. However, superconvergence is, again, lost at the boundary. 
This reduces the total accuracy of the differentiation to 2nd order as with grid x. 
Again, let us note a few explicit examples of differentiation matrices. 

7.2.1 Differentiation Matrix (C, 2pf ) _1 DC^ pt , Grid y 

Again note the unusual, and undesirerable, form of the following matrix. As can be 
seen from Table (7.2.2), again, the superconvergence is lost at the boundary. 


V - 


- 1.99 

3.33 

- 1.60 

.26 

1 

V 

CM 

1 

- le ~ 4 

— 6 e 5 

ift 

1 

1 

le 4 

— 2 e -4 

6 e - 5 

2 e ~ 5 

-.46 

-.14 

.72 

-.13 

- le ~ 3 

— 6 e -4 

— 3 e -4 

— 2 e -4 

7 e " 4 

- le " 3 

3 e -4 

le " 4 

.08 

-.63 

-.02 

.70 

-.12 

- 2 e - 3 

- le " 3 

— 7 e -4 

2 e - 3 

— 4 e -3 

le ' 3 

3 e " 4 

0 

.07 

-.63 

-.02 

.69 

-.12 

— 3 e - 3 

- 2 e - 3 

7 e -3 

-.01 

3 e - 3 

le -3 

0 

0 

.08 

-.63 

-.02 

.69 

-.11 

— 5 e -3 

.02 

-.03 

9 e " 3 

3 e - 3 

0 

0 

0 

.08 

-.63 

-.02 

.69 

-.12 

.05 

-.08 

.02 

7 e -3 

0 

0 

0 

0 

.08 

-.63 

-.02 

.68 

5 e -3 

-.18 

.05 

.02 

0 

0 

0 

0 

0 

.08 

-.63 

-.04 

.94 

-.50 

.12 

.03 

0 

0 

0 

0 

0 

0 

.08 

-.69 

.52 

-.12 

.15 

.07 

0 

0 

0 

0 

0 

0 

0 

-.03 

.43 

- 1.54 

1.00 

.14 

0 

0 

0 

0 

0 

0 

0 

-.20 

2.12 

- 4.25 

2.05 

.29 

0 

0 

0 

0 

0 

0 

0 

.44 

- 4.62 

7.00 

- 2.83 

8 e -3 


It has been mentioned many times that the structure is undesirable. This is due 
to the lack of directional independence of the matrix. One ‘fix’ for this problem is 
given in the next subsection. 

7.2.2 The NPD Differentiation Matrix (C, 2 ^) -1 DC r2p< , Grid y 

NPD denotes No-Preferred-Direction. It is possible to fix one of the problems with 
the above matrices, the problem of directional dependence. That is, a general differ- 
entiation matrix should never treat data coming from the right differently than data 
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coming from the left. This problem can be ‘fixed’ if one simply turns the basis around 
and builds another differentiation matrix. That is, turn all the inner scaling functions 
around and interchange the role of the RHS and LHS boundary functions. The new 
differentiation matrix V can be found from the unreversed differentiation matrix V 
by first flipping the matrix by a middle horizontal axis and then by a middle vertical 
axis and then negating. That is, 


t> = — flipud(fliplr('D )), (55) 

and we build a new NPD differentiation matrix by averaging these two matrices: 

T>npd = + ^)* ( 56 ) 

The following NPD matrix treats data coming from the left in the same manner 
as data coming from the right. Note that this matrix is full. 


- 12.00 

36.96 

- 51.58 

29.30 

- 2.67 

T>npd = 

— 7 e — 4 — 3 e — 4 

- 2 e ” 4 

8 e ” 4 

- lc - 3 

4 e “ 5 

le " 4 

- 4.47 

- 13.12 

29.85 

- 13.47 

1.22 

— 4 e ” 3 

— 2 e ” 3 

- lc ” 3 

4 e ” 3 

— 7 e ” 3 

2 e “ 3 

6 c “ 4 

-.39 

- 9.77 

9.11 

1.62 

-.54 

-.01 

— 6 e ” 3 

- 4 c ” 3 

.01 

-.02 

7 e ” 3 

2 e ” 3 

-.43 

-.43 

- 3.01 

- 3.24 

8.30 

- 1.16 

-.02 

-.01 

.04 

-.07 

.02 

6 e - 3 

-.20 

-.70 

3.44 

- 9.44 

.12 

7.96 

- 1.15 

-.03 

.11 

-.18 

.06 

.02 

-.09 

-.32 

1.09 

.42 

- 7.86 

— 9 e “ 3 

7.94 

- 1.16 

.28 

-.47 

.14 

.04 

-.04 

-.14 

.47 

-.28 

1.16 

- 7.94 

9 c ” 3 

7.86 

-.42 

- 1.09 

.32 

.09 

-.02 

-.06 

.19 

-.11 

.03 

1.15 

- 7.96 

-.12 

9.44 

- 3.44 

.70 

.20 

- 6 e ” 3 

-.02 

.07 

-.04 

.01 

.02 

1.16 

- 8.30 

3.24 

3.01 

.43 

.43 

— 2 c ” 3 

— 7 e " 3 

.02 

-.01 

4 e ” 3 

6 e “ 3 

.01 

.54 

- 1.62 

- 9.11 

9.77 

.39 

1 

<u 

tD 

1 

— 2 e ” 3 

7 e ” 3 

- 4 e ” 3 

le " 3 

2 e ” 3 

4 c ” 3 

- 1.22 

13.47 

- 29.85 

13.11 

4.47 

- le “ 4 

— 4 e ” 4 

lc ” 3 

— 8 e “ 4 

2 e ” 4 

3 e “ 4 

7 e “ 4 

2.67 

- 29.29 

51.58 

- 36.96 

12.00 _ 


Table (7.2.2) gives the accuracy results for two differentiation matrices constructed 
from the well-conditioned quadrature matrix C, 2pt . 

Again, the following table shows that loss of superconvergence at the boundaries 
regardless of the differentiation matrix. Also, note that for grid y that the supercon- 
vergence is not even achieved away from the boundary as it is for grid x. Inner error 
and boundary error are defined in the same manner as for the previous table. 

In conclusion, in this section I have presented a number of differentiation matrices 
along with the orders of accuracy of the differentiation. In all cases the superconver- 
gence which is encountered under periodic boundary conditions is lost. 
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Differentiation 

Grid 

Boundary 

Error 

Inner 

Error 

Total 

Error 

Matrix 

Size 

/i Error 

Ratio 

li Error 

Ratio 

l\ Error 

Ratio 

(C* vt )~ l DC T* 

16 

.0181 


8.70e -4 


.00598 


Grid y 

32 

.00933 

1.94 

2.31e -4 

3.77 

.00162 

3.69 


64 

.00474 

1.96 

6.04e -5 

3.82 

4.26e -4 

3.80 


128 

.00238 

1.99 

1.49e -5 

4.05 

1.089e -4 

3.91 


256 

.00117 

2.03 

3.73e -6 

3.99 

2.72e" 5 

4.00 

NPD of 

16 

.0211 


6.82e~ 4 


.00692 


(Cl pt )~ l DCf pt 

32 

.0116 

1.82 

1.42e -5 

48.0 

.00191 

3.62 

Grid y 

64 

.00606 

1.91 

5.88e -6 

2.41 

5.05e -4 

3.78 


128 

.00310 

1.95 

1.43e -6 

4.11 

1.30e -4 

3.88 


256 

.00156 

1.99 

3.50e -7 

4.09 

3.29e -5 

3.95 


Table 2: Differentiation accuracy for Daubechies wavelets defined on an interval us- 
ing grid points at the boundaries and all other grid points evenly-spaced across the 
interval. 
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8 Other Constructions and Conclusion 

In this paper I have explored the differentiation matrix for Daubechies-type wavelets 
defined on an interval using the boundary construction defined in [3]. The exploration 
has been limited to the wavelet with 2 vanishing moments usually referred to as 
the Daubechies 4 wavelet. It was seen that the superconvergence encountered with 
periodic boundary conditions which was proven in [6] is lost with this boundary 
construction. 

There are other boundary constructions available for Daubechies wavelets, see [1], 
[8], but none of the constructions can maintain the superconvergence. There is a 
very straightforward way to find the maximum accuracy at the boundary for a given 
boundary construction. That is, when the basis functions are orthogonal under trans- 
lation, as they are for Daubechies wavelets, one can simply count the number of basis 
functions which have support overlapping the support of the outermost boundary 
functions. This defines the number of nonzero entries in the first and last rows of the 
middle matrix D referred to here as the derivative projection matrix. All differencing 
in wavelet differentiation is done by D , and, therefore, the differentiation accuracy 
at the boundaries cannot exceed one less than the number of nonzero entries in the 
first and last rows of D. Of course, this argument holds for every row of D, but I am 
particularly concerned with the boundaries here. 

Perhaps it is not possible to construct boundary functions for wavelets which give 
everything. That is, boundary functions which maintain the wavelet structure, which 
maintain the approximation properties across the entire interval and which maintain 
superconvergence. This conjecture is based first on the fact that to date there is not 
a boundary construction which satisfies these criteria for either Daubechies wavelets 
or spline-based wavelets, see [7]. In further support of this conjecture, in a paper 
by Gottlieb, Gustafsson, Olsson, and Strand, [5], it is proven that it is not possible 
to construct boundary functions for a finite element subspace which maintain the 
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superconvergence encountered with periodic boundary conditions if there are charac- 
teristics leaving the domain. Note that the finite element subspace is the same as the 
first order spline wavelet subspace, see [7]. 

In conclusion, under the currently available boundary constructions the supercon- 
vergence, which is encountered under periodic boundary conditions, is not maintained 
at the boundaries. In the case of the D 4 wavelet under the boundary construction 
examined here one gets 4th-order differentiation away from the boundaries and only 
lst-order differentiation at the boundaries producing a scheme with only 2nd-order 
accurate spatial differentiation. It remains to be seen if it is possible to construct 
boundary functions for a wavelet subspace which can maintain superconvergence at 
the boundaries. 
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